$onText
DISE-2024 (Dynamic Integrated Space-Economy) model
File: DISE_2024_03.gms (Social optimum, no-intervention)
Sensitivity analysis: Discount rate
Bongers, A. and Torres, J. L. (2025). Optimal path for orbital debris.
$offText

*	Simulation period
Set t	'Time periods'				/2023*2272/;

Sets
        tfirst(t)   'First period (2023)'
        tlast(t)    'Last period (T)'
    	tnotlast(t) 'All periods but last'
	    ;
        
        tfirst(t)$(ord(t)=1)       =   yes;
        tlast(t)$(ord(t)=card(t))  =   yes;
        tnotlast(t)                =   not tlast(t);

Parameters

* 	Economic parameters

** 	Preferences parameters
	rho	'Intertemporal pure preference rate'		/0.015/
	sigma	'Relative risk-aversion'			/1.5/
    
**  Exogenous growth sources
    	gN0	    'Labor growth rate parameter'			/0.05/
    	Nstar	'Asymptotic population'				/10200/
    	gA0	    'TFP growth rate'				        /0.015/
    	deltaA	'TFP growth rate trend'				/0.001/
    	gQ0	    'Satellite ISTC growth rate'			/0.03/
    	deltaQ	'Satellite ISTC growth rate trend'	/0.005/

**	Initial values economic variables
	k0	    'Initial Earth capital'				/555.6987/         		
    	s0     	'Initial Satellite capital'			/1.1959/  			
    	y0	    'Initial output'				/184.65/
    	x0      'Initial value of satellites destroyed'
    	N0      'Initial population'				/8056/
        
** 	Technological parameters      
	    alpha1 	'Output-capital elasticity'			/0.3479/
    	alpha2 	'Output-satellite elasticity'			/0.0021/ 
    	deltak 	'Capital depreciation rate'			/0.07/
    	deltas 	'Satellite depreciation rate'			/0.15/ 
	    a0	    'Initial value of TFP'				      
	    q0      'Initial value of satellite ISTC'		/1.00/

**  Launch cost
	b0      'Initial proportion of launch cost'     /0.3/
    gb0     'Initial change in launch cost'         /0.0/
    deltab  'Launch cost growth rate trend'         /0.0/

** 	Physical parameters
	deltaf	'Fragments decay rate'				/0.01/
	deltaw  'Dereclit satellites decay rate'        	/0.00015/
	deltaz  'Rocket bodies decay rate'			/0.00015/
	omega 	'Number of fragments per launch'		/4.00/
	xi	    'Fraction of derelict satellites in orbit'	/0.40/
	theta	'Collision risk'				/1.25e-10/
	epsilonw'Fraction of dead satellite breakups'		/0.001/
	epsilonz'Fraction of body rockets breakups'		/0.0012/
	psi	    'Rocket bodies per launch'			/0.60/
	phiw 	'Fragments per derelict satellite breakup'	/44.60/
	phiz    'Fragments per rocket body breakup'		/100.20/
	gammas	'Fragments per operational satellite collision'	/70/
	gammaw  'Fragments per derelict satellite collision'	/70/
	gammaz  'Fragments per rocket body collision'		/70/
	eta	    'Number of satellites per launch'		/13.6/
	v	    'Collision avoidance'				/0/
	mu	    'Mapping parameter'
	
**	Initial values physical variables	
	w0	    'Abandoned derelict satellites 2023'		/3500/
	z0	    'Rocket bodies 2023'				/2050/
	f0	    'Number of fragments larger than 10cm 2023'	/30000/
	d0	    'Total number of debris larger than 10cm 2023'	
	ss0	    'Operational satellites 2023'			/8500/
	xx0	    'Satellites destroyed by collisions 2023'
	ff0	    'Total number of fragments larger than 1cm 2023'
	fff0 	'Total number of fragments larger than 1mm 2023'
	dd0	    'Total number of debris larger than 1cm 2023'
	ddd0 	'Total number of debris larger than 1mm 2023'
	ob1	    'Number of objects larger than 10cm (ESA)'	/30000/
	ob2	    'Number of objects larger than 1cm (ESA)'	/1000000/
	ob3	    'Number of objects larger than 1mm (ESA)'	/130000000/
	nu1	    'Ratio >1cm to >10cm'
	nu2	    'Ratio >1mm to >1cm'
	nu3	    'Ratio >1mm to >10cm'
	;

*	Calculus of key constants and initial values
    	mu	=	ss0/s0;
    	nu1 	=   	ob2/ob1;
    	nu2 	=   	ob3/ob2;
    	nu3 	=   	ob3/ob1;
    	d0 	= 	w0+z0+f0;
	ff0	=	(1+nu1)*f0;
	dd0	= 	w0+z0+ff0;
	fff0	=	nu3*f0;
	ddd0	=	w0+z0+fff0;
	a0	=	y0/(k0**alpha1*s0**alpha2*N0**(1-alpha1-alpha2));
	x0	=	theta*dd0*s0;
	xx0 	=   	theta*dd0*ss0;
	

*   	Timing

Parameters
    	beta(t)	'Discount factor'
    	tval(t)	'numerical value of t';
	    tval(t)	=	ord(t)-1;

Parameters
*   Technology
        A(t)    'TFP'
        Q(t)    'Satellite ISTC'
        gA(t)   'Growth rate of TFG'
        gQ(t)   'Growth rate of Satellite ISTC'
        ;
    
        gA(t)       =   gA0*exp(-deltaA*(tval(t)-1));
        gQ(t)       =   gQ0*exp(-deltaQ*(tval(t)-1));
        a("2023")       =   a0;
        q("2023")       =   q0;
        loop(t,A(t+1)   =   a(t)*exp(gA(t)););
        loop(t,Q(t+1)   =   q(t)*exp(gQ(t)););
Parameters
        b(t)    'Fraction of launch cost'
        gb(t)   'Decline in the launch cost';
        
        gb(t)           =   gb0*exp(-deltab*(tval(t)-1));
        b("2023")       =   b0;
        loop(t,b(t+1)   =   b(t)*exp(gb(t)););
        
    
*   Population
Parameters
    N(t)    'Population';
    N("2023")       =   N0;
    loop(t,N(t+1)   =   N(t)*(Nstar/N(t))**gN0;);

	
* The terminal weight beta(tlast) computation.
	beta(tnotlast(t))	=	power(1+rho,-tval(t));
	beta(tlast(t))		=	(1/rho)*power(1+rho,1-tval(t)-1);
	
Variables

* Economic variables
	SU      'Total summatory discounted utility'
	y(t)	'Output'
   	c(t)    'Consumption'
   	ratcy(t)'Ratio consumption-output'
   	ik(t)	'Investment in Earth capital'
   	is(t)	'Investment in Space capital'
   	h(t)    'Investment in satellite'
   	l(t)    'Launch cost'
   	k(t)    'Earch capital stock'
	s(t)	'Stock of satellites'
	x(t)	'Satellites destroyed by collisions'
	ii(t)   'Investment'
	hh(t)	'Number of new satellites inserted into orbit'
	U(t)	'Instantaneous utility'
	DU(t)   'Discounted utility'
	srk(t)	'Investment rate in Earth capital'
	srs(t)	'Investment rate in Space capital'

* Physical variables
	xx(t)	'Number of destroyed satellites'
	ss(t)   'Number of operational satellites (stock)'
	d(t)	'Number of debris larger than 10cm'
	dd(t)	'Number of debris larger than 1cm'
	ddd(t)	'Number of debris larger than 1mm'
	ll(t)	'Number of launches'
	ds(t)	'Value of derelict satellites per period'
	dss(t)	'Number of derelict satellites per period'
	dso(t)	'Number of derelict satellites abandoned in orbit per period'
	dsb(t)	'Number of derelict satellites breakups'
	dsd(t)	'Number of derelict satellites destroyed by collision'
	w(t)	'Number of derelict satellites abandoned in orbit (stock)'
	rb(t)	'Number of rocket bodies as a function of launches'
	rbb(t)	'Number of rocket bodies breakups'
	rbd(t)	'Number of rocket bodies destroyed by collision'
	z(t)	'Number of rocket bodies (stock)'
	f(t)	'Number of fragments larger than 10cm'
	ff(t)	'Number of fragments larger than 1cm'
	fff(t)	'Number of fragments larger than 1mm'
	;

Equations

* 	Economic variables
	utility(t)      'Instantaneous utility'
	sumdutility	'Summatory of discounted utility'
	production(t)   'Cobb-Douglas production function'
	allocation(t)   'Household choose between consumption and investments'
	totalinvest(t)	'Total investment'
	launch(t)       'Launch cost'
	satel(t)        'Satellite cost'
	accumulak(t) 	'Capital accumulation'
	accumulas(t)	'Satellite accumulation'
	ifinal(t)       'Minimal capital investment in final period'
	hfinal(t)	'Minimal satellite investment in final period'  
	iifinal(t)	'Minimal total investment in final period'
	destroyed(t)	'Value of satellites destroyed by collision'
	ratiocy(t)      'Consumption-output ratio'
	ratioiy(t)	'Investment rate in Earth capital'
	ratiohy(t)	'Investment rate in Space capital'

*	Physical variables
	satellites(t)	'Number of operational satellites'
	deadsat(t)	'Derelict satellites'
	deadsaton(t)	'Derelict satellites abandoned in orbit'
	accumuladso(t)	'Stock of derelict satellites in orbit'
	rocketbodies(t)	'Number of new rocket bodies in orbit'
	accumularb(t)	'Stock of rocket bodies in orbit'
	fragments(t)	'Stock of fragments larger than 10cm in orbit'
	fragments1(t)	'Stock of fragments larger than 1cm in orbit'
	fragments2(t)	'Stock of fragments larger than 1mm in orbit'
	debris(t)	'Stock of debris larger than 10cm'
	debris1(t)	'Stock of debris larger than 1cm'
	debris2(t)	'Stock of debris larger than 1mm'
	newsat(t)	'Number of new satellites inserted into orbit'
	nlaunches(t)	'Number of launches'
	ndeadsat(t)	'Number of derelict satellites per period'
	ndestroyed(t)	'Number of satellites destroyed by collision'
	deadsatb(t)	'Number of derelict satellites breakups'
	deadsatd(t)	'Number of derelict satellites destroyed by collision'
	rbbreakups(t)	'Number of rocket bodies breakups'
	rbdestroyed(t)	'Number of rocket bodies destroyed by collision'
	limit(t)           'Restriction to limit satellites destroyed'
	;


* Equations of the model

** Utility

	utility(t)..	U(t)		=e=	((c(t)/N(t))**(1-sigma)-1)/(1-sigma);
	sumdutility..	SU  	 	=e=	sum(t, beta(t)*U(t)*N(t));

** Production, consumption and investment
	production(t)..    		y(t) 		=e= 	a(t)*(k(t)**alpha1)*(s(t)**alpha2)*(N(t)**(1-alpha1-alpha2));
	allocation(t)..    		c(t) 		=e= 	y(t)-ik(t)-is(t);
	totalinvest(t)..		ii(t)		=e=	    ik(t)+is(t);
	launch(t)..             l(t)        =e=     b(t)*is(t);
	satel(t)..              h(t)        =e=     is(t)-l(t);
	accumulak(tnotlast(t))..  	k(t+1) 		=e=	(1-deltak)*k(t)+ik(t);
	accumulas(tnotlast(t))..	s(t+1)		=e=	(1-deltas)*s(t)+q(t)*(1-b(t))*is(t)-x(t);
	ifinal(tlast)..     		ik(tlast) 	=g= 	(0.02+deltak)*k(tlast);
	hfinal(tlast)..			is(tlast)	=g=	(0.02+deltas)*s(tlast)/q(tlast);
	iifinal(tlast)..		ii(tlast)	=e=	ik(tlast)+is(tlast);
	ratiocy(t)..			ratcy(t)	=e=	c(t)/y(t);
	ratioiy(t)..			srk(t)		=e=	ik(t)/y(t);
	ratiohy(t)..			srs(t)		=e=	is(t)/y(t);

** Damages
	destroyed(t)..			x(t)		=e=	(1-v)*theta*dd(t)*s(t);
	ndestroyed(t)..			xx(t)		=e=	(1-v)*theta*dd(t)*ss(t);

** Satellite sector
	nlaunches(t)..			ll(t)		=e=	mu*h(t)/eta;
	deadsat(t)..			ds(t)		=e=	deltas*s(t);
	ndeadsat(t)..			dss(t)		=e=	deltas*ss(t);
	newsat(t)..			    hh(t)		=e=	eta*ll(t);
	satellites(t)..			ss(t)		=e=	mu*s(t)/q(t);

** Objects in orbit
    	deadsaton(t)..			dso(t)		=e=	xi*deltas*ss(t);
    	accumuladso(tnotlast(t))..	w(t+1)  	=e= 	(1-deltaw)*w(t)-epsilonw*w(t)-theta*(dd(t)+ss(t))*w(t)+dso(t);
    	rocketbodies(t)..		rb(t)       	=e= 	psi*ll(t);
    	accumularb(tnotlast(t))..   	z(t+1)		=e=	(1-deltaz)*z(t)-epsilonz*z(t)-theta*(dd(t)+ss(t))*z(t)+rb(t);
    	fragments(tnotlast(t))..    	f(t+1)   	=e= 	(1-deltaf)*f(t)+omega*ll(t)+phiw*dsb(t)+phiz*rbb(t)+gammas*xx(t)+gammaw*theta*dd(t)*w(t)+gammaz*theta*dd(t)*z(t);
    	fragments1(t)..             	ff(t)    	=e= 	(1+nu1)*f(t);
    	fragments2(t)..             	fff(t)    	=e= 	(1+nu3)*f(t);
    	debris(t)..                 	d(t)		=e= 	w(t)+z(t)+f(t);
    	debris1(t)..                	dd(t)       	=e= 	w(t)+z(t)+ff(t);
    	debris2(t)..                	ddd(t)       	=e= 	w(t)+z(t)+fff(t);
	deadsatb(t)..			dsb(t)		=e=	epsilonw*w(t);
	deadsatd(t)..			dsd(t)		=e=	theta*(dd(t)+ss(t))*w(t);
	rbbreakups(t)..			rbb(t)		=e=	epsilonz*z(t);
	rbdestroyed(t)..		rbd(t)		=e=	theta*(dd(t)+ss(t))*z(t);
    limit(t)..                  x(t)        =l= s(t);

* Bounds
	k.lo(t) = 0.001;
	s.lo(t)	= 0.001;
	c.lo(t) = 0.001;
	y.lo(t) = 0.001;
	is.lo(t)= 0;
	w.lo(t) = 0;
	z.lo(t) = 0;

* Initial conditions
	k.fx(tfirst) 	= 	k0;
	s.fx(tfirst) 	= 	s0;
	y.fx(tfirst)    =   y0;
	d.fx(tfirst) 	= 	d0;
	w.fx(tfirst) 	= 	w0;
	z.fx(tfirst) 	= 	z0;
	f.fx(tfirst)	=	f0;

* Solution options
   	option iterlim = 99900;
   	option reslim  = 99999;
   	option solprint = on;
	option limrow = 0;
	option limcol = 0;
    
model DISE / all /;

solve DISE maximizing SU using nlp;
solve DISE maximizing SU using nlp;


* Export to Excel using GDX utilities
execute_unload "soldise01.gdx"
execute '=gdx2xls soldise01.gdx';



* 3.0% Discount
	rho=0.03;
* The terminal weight beta(tlast) computation.
	beta(tnotlast(t))	=	power(1+rho,-tval(t));
	beta(tlast(t))		=	(1/rho)*power(1+rho,1-tval(t)-1);
	

model DISE2 / all /;

solve DISE2 maximizing SU using nlp;
solve DISE2 maximizing SU using nlp;


* Export to Excel using GDX utilities
execute_unload "soldise102.gdx"
execute '=gdx2xls soldise102.gdx';



* 0.5% Discount
	rho=0.005;
* The terminal weight beta(tlast) computation.
	beta(tnotlast(t))	=	power(1+rho,-tval(t));
	beta(tlast(t))		=	(1/rho)*power(1+rho,1-tval(t)-1);

model DISE3 / all /;

solve DISE3 maximizing SU using nlp;
solve DISE3 maximizing SU using nlp;


* Export to Excel using GDX utilities
execute_unload "soldise103.gdx"
execute '=gdx2xls soldise103.gdx';



* 0.1% Discount
	rho=0.001;
* The terminal weight beta(tlast) computation.
	beta(tnotlast(t))	=	power(1+rho,-tval(t));
	beta(tlast(t))		=	(1/rho)*power(1+rho,1-tval(t)-1);

model DISE4 / all /;

solve DISE4 maximizing SU using nlp;
solve DISE4 maximizing SU using nlp;


* Export to Excel using GDX utilities
execute_unload "soldise104.gdx"
execute '=gdx2xls soldise104.gdx';



